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The potential of Landsat data processing to provide systematic continental scale products has been demon- 
strated by several projects including the NASA Web-enabled Landsat Data (WELD) project. The recent free 
availability of Landsat data increases the need for robust and efficient atmospheric correction algorithms ap- 
plicable to large volume Landsat data sets. This paper compares the accuracy of two Landsat atmospheric cor- 
rection methods: a MODIS-based method and the Landsat Ecosystem Disturbance Adaptive Processing 
System (LEDAPS) method. Both methods are based on the 6SV radiative transfer code but have different at- 
mospheric characterization approaches. The MODIS-based method uses the MODIS Terra derived dynamic 
aerosol type, aerosol optical thickness, and water vapor to atmospherically correct ETM+ acquisitions in 
each coincident orbit. The LEDAPS method uses aerosol characterizations derived independently from each 
Landsat acquisition and assumes a fixed continental aerosol type and uses ancillary water vapor. Validation 
results are presented comparing ETM+ atmospherically corrected data generated using these two methods 
with AERONET corrected ETM+ data for 95 10 kmx 10 km 30 m subsets, a total of nearly 8 million 30 m 
pixels, located across the conterminous United States. The results indicate that the MODIS-based method 
has better accuracy than the LEDAPS method for the ETM+ red and longer wavelength bands. 

© 2012 Elsevier Inc. All rights reserved. 


1. Introduction 

The impact of the atmosphere is variable in space and time and is 
usually considered as requiring correction for quantitative remote 
sensing applications (Liang et al., 2001; Ouaidrari & Vermote, 1999). 
The recent free availability of the United States (U.S.) Landsat data ar- 
chive (Woodcock et al., 2008) has stimulated the development of 
large area Landsat data processing activities. The Web-enabled Land- 
sat Data (WELD) project is systematically generating 30 m composit- 
ed Landsat Enhanced Thematic Mapper Plus (ETM-T) mosaics of the 
conterminous United States (CONUS) and Alaska with the aim of pro- 
viding consistent 30 m data that can be used to derive land cover as 
well as geophysical and biophysical products (Roy et al., 2010). The 
most recent Version 1.5 WELD products include the top of atmo- 
sphere (TOA) reflectance for each of the six reflective wavelength 
Landsat ETM-T bands (Roy et al„ 2011, http://landsat.usgs.gov/ 
WELD.php). The planned Version 2.0 WELD products will be cor- 
rected for atmospheric effects to provide land surface reflectance for 
the approximately 11,000 million and 3100 million 30 m Landsat 
pixels encompassing the CONUS and Alaska respectively. Consistent 
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surface reflectance data are needed to derive geophysical and bio- 
physical products (Vermote et al., 2002). In this paper the accuracy 
of two state-of-the-practice atmospheric correction methods suitable 
for systematic application to large volume Landsat data sets such as 
the WELD data is assessed. 

A number of Landsat atmospheric correction methodologies have 
been developed. Empirical correction methods are mostly variants 
of the dark-object subtraction (DOS) method (Chavez, 1996; Song & 
Woodcock, 2003). In the DOS approach, atmospheric path radiance 
is assumed to be equal to the radiance sensed over dark objects, 
such as dense dark vegetation or deep clear water. After identifica- 
tion, the dark object reflectance is subtracted from the entire image 
for each TOA reflective band. Dark object subtraction based ap- 
proaches do not correct for variations in the atmospheric scattering 
and absorbing constituents across the image or account for multiple 
scattering. Radiative transfer based atmospheric correction ap- 
proaches typically do not have these issues as they model the propa- 
gation of solar electromagnetic radiation through the atmosphere. 
Although radiative transfer based approaches are amenable to sys- 
tematic large volume satellite processing they do require temporally 
and spatially explicit atmospheric characterization data (Vermote et 
al., 2002). 

In this paper, two state-of-the-practice radiative transfer based 
Landsat ETM+ atmospheric correction methods are considered: the 
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Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) 
method (Masek et al., 2006) and a new MODIS-based method described 
herein. Both methods use the 6SV radiative transfer code (available at 
http://6s.ltdri.org/) which has an accuracy better than 1% over a range 
of atmospheric stressing conditions (Kotchenova et al., 2006). The two 
methods differ in the way that they characterize the atmosphere, with 
the greatest difference in the aerosol characterization. The MODIS- 
based method uses MODIS Terra derived aerosol optical thickness 
and dynamic aerosol type to atmospherically correct Landsat ETM+ 
acquisitions in each coincident orbit. The LEDAPS method derives 
the aerosol optical thickness from each Landsat acquisition and inde- 
pendently corrects each acquisition assuming a fixed continental aero- 
sol type. These two correction methods and the data they use are 
described in more detail below. As there are no ground measured surface 
reflectance datasets distributed across the United States with sensor 
footprints and spectral bandpasses similar to Landsat, the accuracies 
of the MODIS-based and LEDAPS atmospherically corrected Landsat 
ETM+ surface reflectance data are quantified by comparison with 
ETM+ surface reflectance derived independently using the 6SV radia- 
tive transfer code parameterized with AERONET sun-photometer re- 
trievals (Dubovik, et al„ 2002; Holben et al., 1998). The AERONET data 
typically enable radiative transfer based atmospheric correction to 2% 
accuracy and have been used previously to validate satellite surface re- 
flectance products (Kotchenova et al., 2006; Vermote et al., 2002). The 
comparisons are undertaken for 10 kmx 10 km Landsat ETM+ spatial 
subsets centered on AERONET sites across the conterminous United 
States, for a year of Landsat ETM+ observations, to provide a compre- 
hensive assessment over a range of land surfaces and atmospheric 
conditions. 


2. Data and pre-processing 

2.1. Landsat ETM+ data 

Level 1 terrain corrected (LIT) ETM+ data were obtained from the 
USGS EROS Landsat Project. The Level 1 T processing includes radiometric 
correction, systematic geometric correction, precision correction using 
ground control chips, and the use of a digital elevation model to correct 
parallax error due to local topographic relief. The CONUS LIT geolocation 
error is less than 30 m (Lee et al., 2004). In this study the six 30 m reflec- 
tive ETM+ wavelength bands were used: blue (0.45-0.52 pm), green 
(0.53-0.61 pm), red (0.63-0.69 pm), near-infrared (0.78-0.90 pm), and 
two middle-infrared bands (1.55-1.75 pm and 2.09-2.35 pm). The six 
bands were converted to top of atmosphere (TOA) reflectance using the 
best available ETM+ calibration coefficients and standard correction for- 
mulae (Grander et al., 2009). A bit-packed band saturation mask was cre- 
ated to define which bands of each pixel were saturated (Roy et al., 2010) 
and two 30 m cloud masks were generated: the Automated Cloud Cover 
Assessment (ACCA) mask (Irish et al., 2006) and a classification tree 
based cloud mask (Roy et al., 2010). 

The 7665 Landsat ETM+ LIT scenes acquired over the contermi- 
nous United States (CONUS) in the period December 1st, 2007 to 
November 30th, 2008 that were used to generate the Version 1 .5 an- 
nual 2008 WELD composite (http://weld.cr.usgs.gov) were consid- 
ered. These data were compared to the geographic locations of the 
119 CONUS AERONET sun-photometer sites. Only Landsat acquisi- 
tions encompassing an AERONET site with reliable cloud-free AERONET 
retrievals on the day of the Landsat 7 overpass were selected. 
This provided 82 Landsat ETM+ acquisitions at 26 AERONET sites 
encompassing surfaces varying from dark vegetation to highly re- 
flective soil (Fig. 1). From these, spatial subsets of 10 kmx 10 km 
centered on the AERONET sites were extracted. A total of 95 
Landsat ETM+ spatial subsets were extracted throughout the year 
(Fig. 2) including dates with snow and with a range of atmospheric 
conditions. 


2.2. Atmospheric characterization data 

The atmospheric characterization data sources for the LEDAPS, 
MODIS-based and AERONET atmospheric corrections of the Landsat 
ETM+ subsets are described below. 

2.2.1. LEDAPS atmospheric characterization data 

The LEDAPS algorithm uses ancillary sea level atmospheric pres- 
sure and water vapor characterization obtained from the National 
Centers for Environmental Prediction (NCEP) and the National Center 
for Atmospheric Research (NCAR), and ozone from the NASA Earth 
Probe Total Ozone Mapping Spectrometer (EP TOMS) (Masek et al., 
2006). The ozone EP TOMS data are defined daily at 1° spatial resolu- 
tion. The sea level atmospheric pressure and water vapor data are de- 
fined every 6 h at 2.5° spatial resolution and their values are 
temporally linearly interpolated to the Landsat overpass time. A static 
global 0.05° digital elevation model was used in the LEDAPS code to 
adjust the atmospheric pressure from sea level to surface level. This 
was achieved by multiplying the sea level pressure with the negative 
exponent of the quotient of the digital elevation and an 8000 m scale 
height (Vermote & Saleous, 2006). 

The aerosol optical thickness (AOT) is retrieved independently 
from each Landsat acquisition using the Kaufman et al. (1997) 
dense dark vegetation (DDV) approach and assuming a fixed “conti- 
nental” aerosol model. LEDAPS DDV pixels are defined as those pixels 
with TOA reflectance p] 0A < 0.15. The mean ETM+ band 1 TOA reflec- 
tance, p] 0A , and the mean band 7 surface reflectance, p s 7 , for DDV 
pixels falling within 1.2 km grids (i.e. 40 30-m Landsat pixels) are 
computed and used to invert for the AOT by using 6S and iteratively 
varying AOT to atmospherically correct p™ until p,«0.33 p 7 . If 
valid DDV targets are not found within a 1.2 km grid cell, the AOT 
value is interpolated by averaging neighboring values within a 7 km 
window. A default AOT (0.06) is used to fill larger gaps. All the atmo- 
spheric characterization parameters are resampled by bilinear inter- 
polation to common 1.2 km grid cells. 

2.2.2. MODIS atmospheric characterization data 

The Landsat 7 ETM+ and MODIS Terra systems are in the same polar 
orbit, with Landsat ETM+ obseivations occurring approximately 
25 min before MODIS Terra nadir observations. The MODIS-based 
method uses the atmospheric characterization data used to correct 
the MODIS Terra TOA reflectance to surface reflectance (Vermote et 
al., 2002). This MODIS-based atmospheric characterization is assumed 
to be the same as for the approximately 25 min earlier Landsat over- 
pass, except for rapidly moving atmospheres. 

The MODIS Terra atmospheric characterization data is defined at 0.05° 
for each MODIS Level 2 granule (approximately 2000 km along track and 
2300 km along scan). The aerosol optical thickness at 550 nm, and the 
aerosol type (low absorption smoke, high absorption smoke, polluted 
urban, and clean urban types) are derived dynamically from the MODIS 
shortwave visible ocean and land bands using an improved non-linear 
version of the Kaufman et al. (1997) dense dark vegetation methodology 
(Vermote & Kotchenova, 2008). The water vapor is derived directly from 
the MODIS near-infrared water vapor bands (typical accuracy 5-10%) 
(Vermote & Kotchenova, 2008), sea-level atmospheric pressure is de- 
fined by NCEP/NCAR 6-hourly Reanalysis data, and NCEP ozone is derived 
from NASA NOAA Total Operational Vertical Sounder (TOVS) ozone re- 
trievals (typical accuracy 0.02 cm-atm). The one arc-second resolution 
ASTER digital elevation model (http://asterweb.jpl.nasa.gov/gdem.asp) 
is used to adjust the atmospheric pressure from sea level to surface 
level in the same manner as the LEDAPS pressure adjustment. 

Fig. 3 shows an example of one day of the 0.05° MODIS Terra de- 
rived aerosol optical thickness (top) and aerosol type (bottom) over 
the conterminous United States (CONUS). A total of three MODIS 
day time overpasses of the CONUS were sensed. The Landsat ETM+ 
sensor has approximately the same nadir ground track as MODIS 
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Fig. 1. Spatial distribution of the 26 AERONET sites across the conterminous United States where there were reliable aerosol measurements at the time Landsat ETM+ overpass 
during the period December 1 2007 to November 30 2008. 


Terra but because of the narrow ETM +15° field of view it only senses 
a swath approximately 185 km across scan. 


2.2.3. AERONET atmospheric characterization data 

The Aerosol Robotic Network (AERONET) is a network of globally 
distributed ground-based sun and sky scanning radiometers that pro- 
vide near-continuous daytime measurements of spectral aerosol opti- 
cal thickness, water vapor, and inversion aerosol products (Dubovik 
et al., 2002; Holben et al., 1998). The AERONET data obtained from 
the AERONET web site (http://aeronet.gsfc.nasa.gov/) include aerosol 
optical thickness in 7-8 narrow spectral bands with center wave- 
lengths from 0.340 pm to 1.640 pm, aerosol volume size distribution 
in 22 bins from 0.05 pm to 15 pm, aerosol complex refractive index 
(real and imaginary components) in four spectral bands with center 
wavelengths of 0.440 pm, 0.657 pm, 0.871 pun and 1.018 pm, the de- 
gree of particle sphericity, and column water vapor (g/cm 2 ). The 
AERONET data include measurement time and date information 
with three data quality levels: Level 1.0 (unscreened), Level 1.5 
(cloud-screened), and Level 2.0 (cloud-screened and quality- 
assured). The Level 2.0 AERONET data are of higher quality but the re- 
trievals are temporally more sparse than the Level 1.5 data. Over the 
119 CONUS AERONET sites, reliable aerosol volume size distribution, 
complex refractive index and the degree of particle sphericity re- 
trievals were not always available at Level 2.0. Consequently, in this 
study, Level 2.0 aerosol optical thickness retrievals and Level 1.5 
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Fig. 2. The monthly number of Landsat ETM+ lOkmxlOkm subsets at the 26 
AERONET sites (Fig. 1) where there were co-incident, reliable aerosol measure- 
ments during the period December 1 2007 to November 30 2008 (an annual total 
of 95 10 km x 10 km subsets). 


aerosol volume size distribution, complex refractive index, and 
degree of particle sphericity retrievals were used. 

The aerosol optical thickness (AOT) was derived as the average of 
the two closest AERONET AOT Level 2.0 retrievals made immediately 
before and after the Landsat overpass time and no more than 30 min 
apart. Similarly, the water vapor was derived as the average of the 
two AERONET water vapor retrievals associated with these two clos- 
est AOT retrievals. The aerosol volume size distribution, complex re- 
fractive index, and degree of particle sphericity retrievals were 
selected as the closet Level 1.5 set occurring the day of Landsat acqui- 
sition with a solar zenith greater than 50° as they are most reliably re- 
trieved when sky radiance measurements are made over a wide range 
of scattering angles (Dubovik, et al., 2002). In order to maximize the 
quality of the AERONET data, only data with a 0.47 pm imaginary 
component of the refractive index less than 0.015 were used. This 
threshold was arbitrarily selected but purposefully quite conserva- 
tive. AERONET retrievals have limitations at low to medium optical 
thickness and CONUS aerosols usually have low absorption (unlike 
for example, African savanna biomass burning aerosols) and imagi- 
nary refractive index values greater than 0.015 are in general suspect 
in those conditions as shown by robustly determined climatological 
aerosol models (Dubovik et al., 2002). 

The surface atmospheric pressure and atmospheric ozone data 
were defined as the average of the LEDAPS and the MODIS values 
for the AERONET site location. The atmospheric pressure and ozone 
data used by the LEDAPS and the MODIS-based methods are essen- 
tially from the same sources, but they are slightly different caused 
by resampling to different spatial resolutions, so an average provides 
an unbiased estimate for the subsequent comparison. 

3. Methods 


3.1. Atmospheric correction 


The 30 m Landsat ETM-I- TOA reflectance for each of the 95 
lOkmxlOkm spatial subsets was atmospherically corrected inde- 
pendently using the AERONET, LEDAPS, and MODIS-based atmo- 
spheric characterizations using 6SV (Kotchenova et al., 2006). The 
atmospheric correction assumes that the surface is Lambertian and 
infinite, and models the TOA reflectance for a given sun-view geome- 
try and spectral band (Kaufman & Sendra, 1988) as: 


p™ 


Patm T 


TJuP S 

1 $atmP 


( 1 ) 


where p TOA is the TOA reflectance, p atm is the atmospheric intrinsic re- 
flectance, T d is the downward atmospheric transmission in the 
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Fig. 3. MODIS Terra derived 0.05° aerosol optical thickness (top) and aerosol type (bottom) characterization data over the conterminous United States for May 17th 2008. On this 
day these data were available for a total of three MODIS day time overpasses of the CONUS. The optical thickness ranged from 0.01 to 1.258. The aerosol types are low absorption 
smoke (blue), high absorption smoke (green), polluted urban (yellow), and clean urban (red). Fill values are shown as black. 


direction of light propagation from the TOA to the surface, T u is the 
upward atmospheric transmission in the direction of light propaga- 
tion from the surface to the sensor, p s is the surface reflectance, and 
s atm is the atmospheric spherical albedo. Note, that the effect of gas- 
eous transmission has been omitted from Eq. (1) for simplicity. 
When atmospheric characterization data are available, p otm , T d , T u , 
and s a[m , for a given sun-view geometry, can be computed using 
6SV and the surface reflectance p s derived. In this study, the geometry 
used was set, as for the LEDAPS code, using the Landsat ETM+ scene 
center solar zenith, the relative azimuth set as the scene center solar 
azimuth, and an arbitrary 3.5° view zenith (half way from nadir to the 
ETM+ swath edge). 

3.1.1. AERONET atmospheric correction 

The 30 m Landsat ETM-I- TOA reflectance was corrected using the 
AEORNET atmospheric characterization data specific to each ETM+ 
subset. The AERONET data were input into the 6SV radiative code 
which was run once per subset to return for each of the six reflective 
bands three coefficients that were used to generate surface reflec- 
tance for each band as: 

. _ p T0A a-b 

(p T0A a—b)c+ 1 

where p s and p™ are the surface and TOA reflectance respectively, 
and the coefficients a, b, c are defined from Eq. (1) as 0 = ^, 
b = and c = s atm . 

The resulting AERONET corrected surface reflectance data are con- 
sidered to provide the surface reflectance “truth”, since the greatest 
uncertainty in atmospheric correction comes from the aerosol 


characterization and the AERONET provides state-of-the-art aerosol 
characterization. 

3.1.2. LEDAPS and MODIS-based atmospheric correction 

The 30 m Landsat ETM+ TOA reflectance was atmospherically 
corrected independently using the LEDAPS and MODIS-based atmo- 
spheric characterization data. This was straightforward for the 
LEDAPS code, which is fully automated and is described in Masek 
et al. (2006), although care was needed to ensure that the saturated 
labeled pixels were read correctly. 

The MODIS-based atmospheric correction required some addi- 
tional development. A multidimensional look-up table was generated 
for each of the four aerosol types by forward modeling with 6SV and 
using the Landsat ETM+ reflective wavelength band characteristics. 
Each look up table returns four atmospheric correction coefficients 
per ETM+ reflective wavelength band. Each look-up table was 
parameterized by ETM+ band number, aerosol optical thickness 
(22 levels), surface atmospheric pressure (7 levels), and sun-view 
geometry (5527 levels covering 0° to 84° sun zeniths, 0° to 69.589° 
view zeniths, and with the azimuthal plane sampled at constant 4° 
scattering angle interval). The ozone, water vapor absorption and 
other gases absorptions were calculated using empirical fits based 
on 6SV forward modeling. 

For each Landsat ETM+ band the four atmospheric correction coef- 
ficients Patm/Td , I'd. T u . and s a(m (defined in Eq. 1) were derived using 
the 0.05° MODIS atmospheric characterization data input into the 
look up table. Fig. 4 (top) illustrates the p olm /T d values in the Landsat 
ETM+ blue band (0.45-0.52 pm) for the conterminous United States 
for May 17th, 2008 computed for the geometiy of a Landsat acquisi- 
tion acquired over path 15 and row 33 near the Chesapeake Bay in 
Maryland. 
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Fig. 4. Atmospheric correction coefficient p am /T d (terms defined in Eq. 1) derived using the 0.05° MODIS atmospheric characterization data for May 17th 2008 computed for the 
geometry of a Landsat acquisition acquired over path 15 and row 33 near the Chesapeake Bay in Maryland (solar zenith 27.2°, view zenith 3.5°, and relative azimuth 129.1°). 


As evident in Figs. 3 and 4 (top) the CONUS MODIS atmospheric 
characterization data have a large amount (typically 50%) of missing 
data, primarily due to clouds at the time of MODIS Terra overpass. 
Consequently, the derived atmospheric correction coefficients 
have a large number of unknown values (black colors in Fig. 4 top). 
Consequently, a spatial interpolation of the four atmospheric correc- 
tion was undertaken. Geostatistical interpolants, such as kriging, are 
computationally expensive (Pringle et al., 2009), spline based inter- 
polants fit to a large surrounding sample data area but the interpo- 
lated values may be outside the range of the sample data, and 
although inverse distance weighting interpolants are computation- 
ally inexpensive they perform poorly for irregular sample data dis- 
tributions (Shepard, 1968). In this study the natural neighbor 
interpolation approach was used as it has elegant properties 
(no tuning parameters, the interpolated values are guaranteed to 
be within the range of the samples used and will pass through the 
input samples and are smooth everywhere except at locations of 
the input samples) (Sibson, 1981 ). Fig. 4 (bottom) shows the natural 
neighbor interpolation of the 0.05° atmospheric correction coeffi- 
cient data illustrated in Fig. 4 (top). The interpolation preserves 
the original data and fills the gaps smoothly. For areas of extensive 
missing data, the interpolation, as with any other method, will be 
less reliable. 

The atmospheric correction coefficients were estimated for each 
30 m ETM+ pixel and band by projecting the 30 m pixel location 
into the natural neighbor interpolated data and bilinear resampling 
the coefficient values from the four neighboring 0.05° values. Since 
the natural neighbor interpolated 0.05° atmospheric correction coef- 
ficients may be less accurate, a quality assessment measure is gener- 
ated as a count of how many (0-4) of the four 0.05° atmospheric 
correction coefficient pixels were interpolated. The MODIS-based 


Landsat ETM+ surface reflectance was then computed, as for the 
AERONET atmospheric correction, as: 

p< = P™/ c i-<* (3) 

P ((p™/c,-c 2 )/c 3 )c 4 + l 1 > 

where p s and p T0A are the surface and TOA reflectance respectively, 
and Ci = Td, c 2 = p a tm/Td, c 3 = T u , and c 4 = s a[m (terms defined for 
Eq. 1). 

3.2. Accuracy assessment methodology 

The LEDAPS and MODIS-based surface reflectance data for each of 
the six ETM+ reflective wavelength bands and for all of the 95 
lOkmxlOkm subsets, were compared pixel-by-pixel to the corre- 
sponding “truth” AERONET surface reflectance data. Only pixels that 
were not saturated in the original LIT data, were not labeled as 
cloudy by either of the ACCA or the classification tree based cloud 
masks, and were not missing due to the Landsat Scan Line Corrector 
issue, were considered. Surface and top of atmosphere (TOA) reflec- 
tance residuals were derived for each pixel as: 

^P i,\ = |P i,\ P i.\.aeronet\ (^) 

a TOA _ I n TOA s I 

^ P i,\ P i,\ P i,\,aeronet\ \*-V 

where p'), A is the LEDAPS or MODIS-based surface reflectance, p T0A i,\ 
is the TOA reflectance, and p s i f \ a eronet is the AERONET surface reflec- 
tance of pixel i for wavelength \, and A p s i: \ is the surface reflectance 
residual for the LEDAPS or the MODIS-based atmospherically 
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corrected reflectance and A p TO \\ is the TOA reflectance residual. The 
residuals are defined in unitless reflectance (scaled 0-1 ). The mean of 
A ,p\ A and of Ap ro '\ Jl was computed for the LEDAPS and MODIS- 
based methods using all the non-saturated, non-cloudy, non- 
missing, pixels in the 95 10 kmx 10 km subsets as: 

Ap\ = (E A P S >>W (6) 

t TOA I A TOA \ 

Ap X = I z_/ A P i,\J/ n \ ( 7 ) 

where A p\ is the mean surface reflectance residual for the LEDAPS or 
the MODIS-based atmospherically corrected reflectance and A p T0 \ is 
the mean TOA reflectance residual and n A is the number of non- 
saturated, non-cloudy pixels for wavelength 

Similarly, the standard deviation of the residuals A p s ,- iA and Ap m \\ 
were computed as: 

ct aPa = (Vi,x-Ap\)/n A j (8) 

Oa p rm K = ^ (Ap™ u -Ap™ x )/n A j (9) 

where a Ap s A and 0 ^™^ are the standard deviation of the residuals for 
the LEDAPS (or MODIS-based) surface reflectance and the TOA reflec- 
tance respectively and the other terms are defined earlier. 

In order to be able to inter-compare the residuals between spec- 
tral bands, mean reflectance normalized residuals were derived as: 


A*p\ = « 


A P s , 


P A ,aeronet 


a n TOA 

A*p T0 \ = — 

P A ,aeronet 


( 10 ) 


(ii) 


where A *p s A and A *p ro \ are the mean reflectance normalized resid- 
ual for the LEDAPS (or MODIS-based) surface reflectance and the TOA 
reflectance respectively, A p\ and A p TO \ are defined in Eqs. (6) and 
(7) respectively, and p s is the mean of the AERONET corrected 
reflectance at wavelength A using all the non-saturated, non-cloudy, 
non-missing, pixels in the 95 10 kmx 10 km subsets. 

Reflectance scatter plots and simple linear regressions between 
the TOA (y axis) and AERONET surface reflectance (x axis) were gen- 
erated from all the non-saturated, non-cloudy, non-missing, pixels in 
the 95 10 kmx 10 km subsets to assess the impact of the atmosphere 
on each ETM+ reflective wavelength band. Then reflectance scatter 
plots and simple linear regressions between the LEDAPS or MODIS- 
based surface reflectance (y axis) plotted against the AERONET sur- 
face reflectance (x axis) were used to assess the variability and the 
bias in the corrected reflectance. 


4. Results 


A total of 7,605,732, 7,605,785, 7,605,338, 7,606,184, 7,605,977, and 
7,606,020 30 m pixels in Landsat ETM+ bands 1 (0.45-0.52 pm), 2 
(0.53-0.61 pm), 3 (0.63-0.69 pm), 4 (0.78-0.90 pm), 5 (1.55-1.75 pm), 
and 7 (2.09-2.35 pm) respectively that were not saturated in the original 
LIT data and that were not labeled as cloudy by either of the ACCA or the 
classification tree based cloud masks were extracted from the 95 
10 kmx 10 km subsets. These numbers differ between bands because 
the Landsat ETM+ saturation varies spectrally (Markham et al, 2006). 


Figs. 5 and 6 show reflectance scatter plots of TOA reflectance ver- 
sus AERONET surface reflectance (left column), LEDAPS surface re- 
flectance versus AERONET surface reflectance (middle column), and 
MODIS-based surface reflectance versus AERONET surface reflectance 
(right column). The solid lines show ordinary least squares linear re- 
gression fits of these data. If the LEDAPS or MODIS surface reflectance 
data were corrected perfectly then they would have equal value as 
the AERONET surface reflectance data and all points in the middle 
and right columns of Figs. 5 and 6 would fall on the dotted 1:1 lines. 
For the LEDAPS and MODIS-based surface reflectance plots (middle 
and right columns) the regression lines are constrained to pass 
through the origin to examine which method provides generally clos- 
est agreement with the AERONET surface reflectance (i.e. slopes clos- 
er to unity). The frequency of occurrence of the nearly 8 million 
reflectance values is shown with a rainbow color scale (red most fre- 
quent, purple least frequent). The regression coefficients and good- 
ness of fit (R 2 ) values are shown on the figures — all of the 
regressions and R 2 values were significant at the 99% confidence level. 

Fig. 5 shows plots for the three Landsat ETM+ visible wavelength 
bands: band 1 (blue, 0.45-0.52 pm) top row, band 2 (green, 
0.53-0.61 pm) middle row, and band 3 (red 0.63-0.69 pm) bottom 
row. These are the shortest wavelength bands and atmospheric scat- 
tering is expected to be greatest at shorter wavelength. This is evident 
in that the slopes of the regression lines for the blue, green and red 
TOA reflectance against the AERONET surface reflectance are 0.74, 
0.79 and 0.86 respectively (Fig. 5, first column). These three regres- 
sion slopes are less than unity because Rayleigh and aerosol backscat- 
ter into the sensor adds to the TOA signal at low reflectance ranges 
and aerosol absorption attenuates the TOA signal at higher reflec- 
tance. The effect of the LEDAPS and the MODIS-based atmospheric 
corrections is to provide surface reflectance that is closer to the AERONET 
surface reflectance (Fig. 5, middle and right columns), with both 
methods providing surface reflectance estimates that overall slightly un- 
derestimate the AERONET surface reflectance (regression slopes vaiying 
from 0.947 to 0.993). The MODIS-based surface reflectance shows slight- 
ly better or comparable linear relationships with the AERONET surface 
reflectance (higher R 2 values and slopes closer to unity) than the LEDAPS 
surface reflectance in the green and red bands, but a slightly worse linear 
relationship in the blue band. 

Fig. 6 shows plots for the longer wavelength Landsat ETM+ 
bands: band 4 (NIR, 0.78-0.90 pm) top row, band 5 (middle-lR, 
1.55-1.75 pm) middle row, and band 7 (middle-IR, 2.09-2.35 pm) 
bottom row. In this longer-wavelength range, the effect of the at- 
mosphere is less apparent and the TOA reflectance is about 0.92, 
0.95 and 0.89 of the surface reflectance respectively (Fig. 6, first 
column). Consequently the difference between the MODIS-based 
and LEDAPS surface reflectance is less apparent although for all 
three bands the MODIS-based surface reflectance shows slightly 
better linear relationships with the AERONET surface reflectance 
(slopes closer to unity) than the LEDAPS surface reflectance 
(Fig. 6, middle and right columns). Both methods provide surface 
reflectance estimates that slightly overestimate the AERONET sur- 
face reflectance, particularly the LEDAPS results (regression slopes 
varying from 1.036 to 1.054), compared to the MODIS results 
(regression slopes vaiying from 1.007 to 1.014). It is unclear 
what the causes of these biases are. 

Table 1 summarizes, for each Landsat ETM+ reflective wavelength 
band, the mean LEDAPS and MODIS-based surface reflectance resid- 
uals (Eq. 6) and also for reference the mean TOA reflectance residuals 
(Eq. 7). The standard deviations of these residuals (Eqs. 8 and 9) are 
summarized in parentheses. The mean TOA reflectance residuals 
show the significantly increasing impact of the atmosphere with decreas- 
ing wavelength. For example, in the ETM+ blue band (0.45-0.52 pm) the 
mean TOA reflectance residual is 0.0669, about two orders of magnitude 
greater than the mean LEDAPS or mean MODIS-based blue band surface 
reflectance residuals. 
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Fig. 5. Visible wavelength band Landsat ETM+ reflectance scatter plots of TOA reflectance versus AERONET surface reflectance (left column), LEDAPS surface reflectance versus 
AERONET surface reflectance (middle column), and MODlS-based surface reflectance versus AERONET surface reflectance (right column). The solid lines show ordinary least 
squares linear regression fits of these data. For the LEDAPS and MODlS-based surface reflectance plots (middle and right columns) the regression lines are constrained to pass 
through the origin. The dotted lines are 1 : 1 lines superimposed for reference. The points show the 30 m pixel values colored with a rainbow scale to illustrate the frequency of 
the reflectance values (red most frequent, purple least frequent). 


The most striking result apparent in Table 1 is that the mean and 
standard deviation of the MODlS-based surface reflectance residuals 
are smaller than those of the LEDAPS surface reflectance residuals in 
all the Landsat ETM+ bands except for the green and blue bands. 
The green band mean residuals are the same for the LEDAPS and 
MODlS-based methods but the MODIS residuals have smaller stan- 
dard deviation. The spectral variation in the LEDAPS and MODlS- 
based surface reflectance residuals is complex and is driven by the 
spectral variation in atmospheric contamination (which generally de- 
creases with wavelength) and also by the spectral properties of the 
surface, whereby for example, healthy vegetation has low red reflec- 
tance and high near-infrared reflectance. To investigate this, the 
mean reflectance normalized residuals (Eqs. 10 and 11), i.e„ the 
mean residuals expressed as percentages of the mean AERONET sur- 
face reflectance, were computed and are summarized in Table 2, 

In the blue band the LEDAPS and MODlS-based mean reflectance 
normalized residuals are 11.8% and 13.5% respectively (Table 2). 
This suggests that both the LEDAPS and MODlS-based methods pro- 
vide unreliable blue band atmospheric correction, which is a known 
atmospheric correction issue due to the high atmospheric sensitivity 
at this wavelength (Vermote & Kotchenova, 2008) and is evident in 
that the mean reflectance normalized TOA residual is nearly 150% 
(Table 2, right column). Both methods provide green band mean re- 
flectance normalized residuals of 5.7% which is nearly seven times 
less than the mean green band reflectance normalized TOA residual 
(39.9%). In the red band the LEDAPS mean reflectance normalized re- 
sidual (5.9%) is greater than that of MODIS (4.2%) and both mean 


residuals are about four times smaller than the mean TOA residual 
(19.3%). In the longer wavelength bands the LEDAPS mean reflec- 
tance normalized residuals are more than two to more than three 
times greater than the MODlS-based mean reflectance normalized 
residuals. 

The results summarized in Tables 1 and 2 indicate that the MODlS- 
based atmospheric correction approach is generally more accurate 
and robust than the LEDAPS approach. For all the Landsat ETM+ 
bands (Table 2), the TOA mean reflectance normalized residuals are 
significantly greater than the MODIS mean reflectance normalized re- 
siduals, illustrating the need for Landsat atmospheric correction. Sim- 
ilarly, the TOA mean reflectance normalized residuals are greater than 
the LEDAPS mean reflectance normalized residuals, except for the 
near-infrared Landsat band 4 which has a 4.8% LEDAPS mean reflec- 
tance normalized residual compared to the 4.1% TOA band 4 mean re- 
flectance normalized residual. This does not mean that atmospheric 
correction is generally not needed for Landsat band 4, and we note 
that band 4 is also the band with the smallest relative difference be- 
tween MODIS mean reflectance normalized residuals compared to 
the TOA mean reflectance normalized residual. It is well established 
that near-infrared wavelengths are particularly susceptible to water 
vapor contamination (Vermote & Kotchenova, 2008; Vermote & 
Saleous, 2006) and so it is likely that this LEDAPS band 4 discrepancy 
may be due to spatial and/or temporal resolution mismatches, for cer- 
tain Landsat acquisitions, between the atmospheric water vapor con- 
tent at the time of Landsat overpass and the 2.5° six hour NCEP/NCAR 
water vapor data used by the LEDAPS correction approach. The 
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Fig. 6. Infrared reflective wavelength band Landsat ETM+ reflectance scatter plots (as Fig. 5). 


MODIS approach uses 0.05° water vapor retrieved from the MODIS 
overpass which is not only of higher spatial resolution but also nearly 
coincident with the Landsat acquisition time. In addition, however, 
the way that the different atmospheric characterization data are 
interpolated to match Landsat acquisition times and locations is also 
likely to be influential. This is demonstrated in Fig. 7 which shows a 
spectral plot of the mean MODIS-based surface reflectance residuals 
(Eq. 6) considering all the 30 m pixels (open symbols, same results 
as Table 1 ) and considering only 30 m pixels where no 0.05° correc- 
tion coefficients were natural neighbor interpolated (filled symbols). 
For the almost 8 million 30 m pixels used in this analysis, approxi- 
mately 38% had at least one natural neighbor interpolated 0.05° cor- 
rection coefficient pixel. As expected, the mean MODIS-based 
surface reflectance residuals illustrated in Fig. 7 are consistently 
higher, by about 5-9% for ETM+ bands 1-4, and similar in ETM+ 
band 5, when no natural neighbor interpolated values are used. 


Conversely, the mean MODIS-based surface reflectance residuals for 
ETM+ band 7 are higher by 8% when no natural neighbor interpolat- 
ed values are used. The results in Fig. 7 suggest the utility of providing 
a quality assessment band with the MODIS-based atmospheric correc- 
tion results that count, for example, how many of the four 0.05° 
MODIS atmospheric correction coefficient pixels were bilinearly inter- 
polated in the generation of a 30 m ETM+ surface reflectance value. 

It is beyond the scope of this paper to analyze the impact of the at- 
mospheric correction methodologies on higher level products that 
may be derived from Landsat surface reflectance. However, the Nor- 
malized Difference Vegetation Index (NDVI), defined as the near- 
infrared minus the red reflectance divided by their sum, is widely 
used and, for this reason, is included as a WELD product band (Roy 
et al., 2010). To examine the impact on NDVI of the two atmospheric 
correction approaches, NDVI residuals and mean residuals were com- 
puted as Eqs. (4) and (6) respectively (but using the NDVI instead of 


Table 1 

The mean residuals of the LEDAPS surface reflectance and the MODIS-based surface re- 
flectance (Eq. 6), and the mean residuals of the TOA reflectance (Eq. 7), for each Land- 
sat ETM+ reflective wavelength band. The standard deviations of the residuals are 
shown in parenthesis. 


Band 

LEDAPS surface 
reflectance A ff\ and 

O'Ap'A 

MODIS-based surface 
reflectance A and Oa^a 

TOA reflectance 
A p™A and CTip>™ A 

1 

0.0053 (0.0058) 

0.0060 (0.0052) 

0.0669 (0.0111) 

2 

0.0039 (0.0043) 

0.0039 (0.0034) 

0.0276 (0.0086) 

3 

0.0042 (0.0039) 

0.0030 (0.0026) 

0.0138 (0.0072) 

4 

0.0100 (0.0079) 

0.0041 (0.0039) 

0.0085 (0.0063) 

5 

0.0056 (0.0049) 

0.0015 (0.0015) 

0.0068 (0.0040) 

7 

0.0051 (0.0051) 

0.0016 (0.0016) 

0.0105 (0.0079) 


Table 2 

The mean AERONET surface reflectance, and the mean reflectance normalized residuals 
of the LEDAPS surface reflectance and the MODIS-based surface reflectance (Eq. 10), 
and the mean reflectance normalized residuals of the TOA reflectance (Eq. 1 1 ). 


Band 

Mean AERONET 
surface reflectance 

(p S \,aeronet) 

Mean 
reflectance 
normalized 
residual for 
LEDAPS ?(A *p\) 

Mean reflectance 
normalized 
residual for 
MODIS-based 
(A*P s a) 

Mean 
reflectance 
normalized 
residual for 
TOA (A "p TO \) 

1 

0.0447 

11.8% 

13.5% 

149.7% 

2 

0.0692 

5.7% 

5.7% 

39.9% 

3 

0.0714 

5.9% 

4.2% 

19.3% 

4 

0.2092 

4.8% 

2.0% 

4.1% 

5 

0.1550 

3.6% 

1.0% 

4.4% 

7 

0.0986 

5.2% 

1.6% 

10.6% 
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Fig. 7. Spectral plot of the mean MODIS-based surface reflectance residuals (Eq. 6) con- 
sidering all ~8 million 30 m pixels (open symbols, same results as Table 1 ) and consid- 
ering only the ~5 million 30 m pixels where no 0.05° correction coefficients were 
natural neighbor interpolated (filled symbols). 

spectral reflectance in these formulae), and were then normalized 
with respect to the mean AERONET NDV1 as Eq. (10), using NDVI de- 
rived from the AERONET, LEDAPS and MODIS-based bands 3 and 4 
surface reflectance data. The mean NDVI normalized residuals were 
3.1% and 6.3% for the MODIS-based method and the LEDAPS method 
respectively, i.e., similar to the maximum of the normalized bands 3 
and 4 surface reflectance residuals for each method reported in 
Table 2. For reference, the mean TOA NDVI normalized residual, com- 
puted as Eqs. (5), (7) and (11) using NDVI derived from TOA surface 
bands 3 and 4 reflectance, was 16.8%. Evidently the impact of atmo- 
spheric correction on NDVI is quite important which has been ob- 
served by many other researchers (Miura et al., 2001; Slater & 
Jackson, 1982; Teillet et al„ 1997). 

5. Discussion and conclusions 

This study compared the atmospheric correction accuracy of the 
Landsat Ecosystem Disturbance Adaptive Processing System 
(LEDAPS) and a new MODIS-based atmospheric correction method. 
A comprehensive validation was undertaken. A total of nearly 8 mil- 
lion Landsat ETM+ pixels were atmospherically corrected using the 
two methods and compared with AERONET corrected equivalents. 
These data were extracted across the conterminous United States 
from 82 Landsat ETM+ acquisitions sensed from December 2007 to 
November 2008 at 26 AERONET sites. The results indicate that the 
MODIS-based method has overall higher accuracy than the LEDAPS 
method for all the ETM+ bands except the green band, where the re- 
sults for the two methods are comparable, and the blue band where 
both the LEDAPS and MODIS-based atmospheric correction methods 
performed less reliably. In the longer wavelength reflective bands 
the LEDAPS atmospheric correction method performed less reliably 
than the MODIS-based method by a factor of about two to three. 
The accuracy level provided by the MODIS-based Landsat ETM+ at- 
mospheric correction method is comparable to that provided by the 
operational MODIS atmospheric correction algorithm implemented 
over a diverse range of surfaces and atmospheres (Vermote & 
Kotchenova, 2008; http://modis-sr.ltdri.org/). This degree of agree- 
ment is not surprising because, although MODIS has superior spectral 
and radiometric characteristics compared to the Landsat ETM+, at 
scan edge the MODIS atmospheric path length is nearly double the 
ETM+ path length due to the difference in the field of view of the 
MODIS (110°) and Landsat ETM+ (15°) instruments. 

The LEDAPS and the MODIS-based methods both use the 6SV radi- 
ative transfer code, but differ in the way that they characterize the 


atmosphere. The MODIS-based method uses MODIS Terra derived 
aerosol optical thickness, aerosol type and water vapor to atmospher- 
ically correct Landsat ETM+ acquisitions in each coincident orbit. The 
LEDAPS method derives the aerosol optical thickness from each Land- 
sat acquisition and independently corrects each image assuming a 
fixed continental aerosol type, and uses 2.5° spatial resolution water 
vapor defined eveiy 6 h from the National Centers for Environmental 
Prediction (NCEP) and the National Center for Atmospheric Research 
(NCAR). The sensitivity of 6S to different atmospheric parameteriza- 
tions is non-linear and difficult to extrapolate to all combinations of 
characterizations (Vermote & Saleous, 2006). The MODIS Terra in- 
strument senses a much larger swath (-2300 km) than the ETM+ 
(-185 km) which provides more opportunities for dense dark vegeta- 
tion (DDV) target identification and the improved version of the 
MODIS DDV approach allows for brighter vegetation targets for aero- 
sol retrieval. Further, MODIS has superior spectral and radiometric 
characteristics compared to the Landsat ETM+ instrument and the 
MODIS near-infrared bands can be used for water vapor retrieval 
(Vermote & Kotchenova, 2008). For these reasons, the MODIS-based 
method should provide more reliable atmospheric characterization 
than the LEDAPS approach. However, the MODIS spatial resolution 
(500 m) is much coarser than the ETM+ spatial resolution (30 m) 
at the wavelengths used for DDV target identification, and the 
MODIS atmospheric characterization describes the atmosphere ap- 
proximately 25 min after the Landsat ETM+ overpass. Consequently, 
for example, small and spatially fragmented DDV targets less than 
500 m in dimension and dynamic aerosols may be better defined 
from the ETM+ acquisition itself under the LEDAPS approach. In ad- 
dition, the way that the different atmospheric characterization data 
are interpolated to match Landsat acquisition times and locations is 
likely to influence the atmospheric correction accuracy. Sensitivity 
analyses are required to investigate the results reported in this 
paper in more detail, the multiple influencing factors behind them, 
and the impacts of error in surface reflectance on higher level de- 
rived Landsat products and on applications such as land cover 
characterization. Certainly however, improved atmospheric char- 
acterization at the time of satellite overpass will result in more re- 
liable radiative transfer based atmospheric corrections. For this 
reason further work to investigate a fusion of Landsat image- 
based aerosol retrievals and MODIS-based atmospheric characteri- 
zation data is also recommended. 

Finally, we note that although the MODIS-based atmospheric cor- 
rection has generally better performance than the image-based 
LEDAPS approach, the LEDAPS approach can be applied to the historic 
Landsat Thematic Mapper archive (available since 1982) and also to 
future Landsat sensors regardless of the availability of the next gener- 
ation of operational moderate spatial resolution global polar-orbiting 
remote sensing systems (Roy et al., 2008; Wulder et al„ 2011). 
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